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Abstract 

A thermodynamically self-consistent Ornstein-Zernike approximation (SCOZA) 
is applied to a fluid of spherical particles with a pair potential given by a hard-core 
repulsion and a Yukawa attractive tail w(r) = — exp[— z(r — l)]/r. This potential 
allows one to take advantage of the known analytical properties of the solution to 
the Ornstein-Zernike equation for the case in which the direct correlation function 
outside the repulsive core is given by a linear combination of two Yukawa tails 
and the radial distribution function g(r) satisfies the exact core condition g(r) = 
for r < 1. The predictions for the thermodynamics, the critical point, and the 
coexistence curve are compared here to other theories and to simulation results. 
In order to unambiguously assess the ability of the SCOZA to locate the critical 
point and the phase boundary of the system, a new set of simulations has also 
been performed. The method adopted combines Monte Carlo and finite-size scaling 
techniques and is especially adapted to deal with critical fluctuations and phase 
separation. It is found that the version of the SCOZA considered here provides 
very good overall thermodynamics and a remarkably accurate critical point and 
coexistence curve. For the interaction range considered here, given by z = 1.8, the 
critical density and temperature predicted by the theory agree with the simulation 
results to about 0.6%. 
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1 Introduction 



After applying their version of thermodynamic perturbation theory to square-well and 
Lennard- Jones fluids, John Barker and Doug Henderson characterized it as a "successful 
theory of liquids" ][]]. And so it was. When tested against simulation results it proved to 
be impressively accurate at liquid-state densities and temperatures, unlike some versions 
of thermodynamic perturbation theory that had preceded it. And it bypassed the trou- 
bling lack of thermodynamic self-consistency associated with the direct use of the radial 
distribution functions obtained from the integral-equation theories then available, as well 
as yielding thermodynamic results as good or better than the best results obtainable from 
such integral equations. 

These positive features became hallmarks of successful thermodynamic perturbation 
theories for simple fluids and were shared by the versions @ that followed the Barker 
and Henderson work as well as an alternative perturbative approach set forth somewhat 
earlier by Hauge and Hemmer |§ that was based on using the inverse range of the attrac- 
tive interaction rather than its strength as a perturbation parameter. Integral-equation 
approaches with improved self-consistency were also developed subsequently to yield ac- 
curate liquid-state thermodynamics [f|. 

Unfortunately, the accuracy of all these approaches begins to decrease substantially 
as one leaves the liquid-state region located slightly above the triple point in temperature 
and follows the liquid-gas coexistence curve in the density-temperature plane up to the 
critical region. In particular, the shape of the coexistence curve and location of the critical 
point are not accurately reproduced, nor are related critical parameters. In the case of 
the perturbation theories, it is not hard to understand why this is so. All of them are 
mean-field-like in nature, associated with coexistence curves that are quadratic close to 
the critical point, whereas the true coexistence curve is very nearly cubic. That is, in 
these theories one finds near the critical point a coexistence curve of the form 

T c -T*A\p-p c \ x , x = 2, (1) 

where p c and T c are the critical values of number density p and absolute temperature 
T, and A is a constant. In contrast, in an exact treatment, one would expect to find x 
close to 3. In these theories the resulting T c is usually more than 5% too high and the 
critical compressibility factor (P/ pksT) c is usually more than 10% too high. Here P is 
the pressure, and ks is the Boltzmann constant. 

The thermodynamics associated with the radial distribution function g(r) obtained 
form various integral-equation approaches cannot be so neatly categorized. However, in 
the cases in which there are substantial discrepancies between the several paths available 
for obtaining thermodynamics from g(r), the most reliable and accurate coexistence be- 
havior is often obtained from evaluating the thermodynamics through the excess internal 
energy expressed in terms of an integral over the pair potential w(r) weighted by g(r). For 
continuum-fluid models the resulting critical behavior is typically mean-field like in the 
cases that we have studied, and thus subject to the same deficiencies as one approaches 
the critical region. In some integral-equation approaches that have been developed in or- 
der to insure a certain degree of thermodynamic consistency, the description of the critical 
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region and of the phase diagram appears to be more problematic: for instance, the mod- 
ified hypernetted chain (MHNC) theory || is indeed able to predict quite satisfactorily 
the liquid and the vapor branches of the coexistence curve of a simple fluid at low enough 
temperature, but it fails to converge close to the critical point, so that the two branches 
remain unconnected, and the position of the critical point is not given directly by the 
theory, but must be determined by extrapolation || [7| . The same kind of behavior || |[ 
is found also for the HMSA integral equation [the acronym coming from the fact that the 
theory || interpolates between the hypernetted chain (HNC) and the soft mean spherical 
approximation (SMSA)]. 

The self-consistent Ornstein-Zernike approximation (SCOZA) we consider here is not 
mean-field-like, and it remains highly accurate as one goes from liquid-state conditions 
to critical-point conditions. In particular the power x in Eq. ([[]) was recently shown 
analytically to be given in the SCOZA by exactly 20/7 ||10|| . And as we discuss in this 
paper, in the hard-core Yukawa fluid (HCYF) T c appears to be within 0.6% of its value 
as estimated by our simulation results. (Similarly, in recent three-dimensional lattice-gas 
studies JO], [12[| the SCOZA T c was found to be within 0.2% of its estimated exact value). 
As described elsewhere [pi], [12|] the scaling behavior of the SCOZA thermodynamics is 
somewhat different from the simple scaling one expects to see in the exact thermodynam- 
ics, although those differences only begin to appear clearly when p and T are within less 
than 1% of their critical values. Closer to the critical point, the effective exponents defined 
above T c approach spherical-model values as the critical point is approached, whereas the 
exponents defined below T c do not. The exponents are discussed in Sec. 3. 

The SCOZA was proposed some time ago by H0ye and Stell [13|, [L4j but fast and 
accurate algorithms for evaluating its thermodynamic predictions were developed only 
recently fll , |12| 15 1. A sharp assessment of its accuracy for the HCYF could not be 



made on the basis of existing simulations, and for that reason our study here includes 
new Monte Carlo (MC) results exploiting finite-size scaling (FSS) techniques (16| . 

We have chosen the HCYF pair potential as the first of the continuum-fluid potentials 
to be considered in our studies of the SCOZA for several reasons. First, it embodies the 
two key features one requires in an off-lattice pair potential in order to consider both 
the liquid state and liquid-gas criticality — a highly repulsive core and an attractive well. 
Second, the HCYF proves to be particularly convenient to analyze using the SCOZA 
(the square- well fluid is far less convenient in this regard). Third, the functional form 
of the hard-core Yukawa potential makes it appropriate as a generic solvent-averaged 
interaction potential between polyelectrolytes and colloids as well as a generic simple- 
fluid pair potential. For this reason it seems particularly useful to have an accurate 
theory for both the structure and thermodynamics of the HCYF, which has already been 
the subject of a number of previous studies. We shall make contact with several of those 
here. 

The paper is organized as follows: in Sec. 2 we describe the theory and present some 
details of the method for the system under study, in Sec. 3 our results are shown and 
a comparison with other theories and simulation results is made, and in Sec. 4 our con- 
clusions are drawn. The treatment of the hard-sphere gas and the main features of the 
MC-FSS simulation method are summarized respectively in Appendix A and Appendix B. 
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2 Theory 




Here we consider a fluid of spherical particles interacting via a two-body potential v(r) 
which is the sum of a singular repulsive hard-sphere contribution and an attractive tail 
w(r) < 0. The expression for v(r) is then 

r < 1 

(2) 

r > 1 , 

where the hard-sphere diameter has been set equal equal to unity. As is customary 
in integral equation theories of fluids, the present approach introduces an approximate 
closure relation for the direct correlation function c(r) which, once supplemented with the 
exact Ornstein-Zernike equation involving c(r) and the radial distribution function g(r), 
yields a closed theory for the thermodynamics and the correlations of the system under 
study. The basic requirement we want to incorporate in the SCOZA is the consistency 
between the compressibility and internal energy route to the thermodynamics. According 
to the compressibility route, the thermodynamics stems from the reduced compressibility 
Xred as determined by the sum rule 

Xrcd = 1 - pc(k = 0) ' (3) 

where c(k) is the Fourier transform of the direct correlation function and p is the number 
density of the system. In the internal energy route the key to the thermodynamics is 
instead provided by the excess internal energy as given by the integral of the interaction 
weighted by the radial distribution function: 

/+oo 
drr 2 w(r)g(r) , (4) 

where u is the excess internal energy per unit volume and we have taken into account 
that g(r) vanishes for r < 1 due to the hard-core repulsion. In the following we will refer 
to the "excess internal energy" simply as the "internal energy" . If Xrcd and u come from 
a unique Helmholtz free energy it is straightforward to find that one must have 

where j3 = l/(/csT), T being the absolute temperature, and ks the Boltzmann constant. 
While this relation is of course satisfied by the exact compressibility and internal energy, 
this is not the case with those predicted by most integral equation theories. In order to 
cope with this lack of thermodynamic consistency, we consider the following closure to 
the Ornstein-Zernike equation: 

' g(r) = r < 1 , 

(6) 

k c(r) = c ns {r) + K(p, (3)w(r) r > 1 , 
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where CHs( r ) is the direct correlation function of the hard-sphere fluid, and K(p,/3) is a 
function of the thermodynamic state of the system. In Eq. @ the approximation clearly 
lies in the simple form of c(r) outside the repulsive core. The closure above resembles 
the one used in the approximation known as both the lowest-order gamma-ordered ap- 
proximation (LOG A) [[17] and the optimized random phase approximation (ORPA) fll8[ . 
However, while in the LOGA/ORPA one has K(p,/3) = -p, in Eq. (§) K(p,/3) is not 
fixed a priori, but instead must be determined so that the thermodynamic consistency 
condition ([5]) is satisfied. This gives rise to a partial differential equation (PDE) for the 
function K(p,/3), provided an expression for the hard-sphere direct correlation function 
c Hs( r ) is given. The most popular parameterization for CHs( r ) i n the fluid region is due 
to Verlet and Weis [fL9[j . Another choice that yields the same thermodynamics as Verlet- 
Weis, and that we find convenient in view of the calculations performed in this work, is 
originally due to Waisman [20]. It was subsequently extended analytically by H0ye and 
Stell |]2TJ and explored in some detail by Henderson and coworkers ||22||. It amounts to 



assuming that the function CHs( r ) outside the repulsive core has a one- Yukawa form, so 
that for the hard-sphere system we have: 



' 9us(r) = 



r < 1 



exp[-zi(r - 1)] 
c H s(^)=^i r>l. 



(7) 



The Ornstein-Zernike equation supplemented by Eq. fl7|) can be solved analytically in 
terms of the amplitude K\ and the inverse range z\ of CHs( r )- These can be in turn 
determined as a function of the density by requiring, as in the Verlet- Weis parameteriza- 
tion, that both the compressibility and the virial route to the thermodynamics give the 
Carnahan-Starling equation of state. The basic features of the calculation are recalled in 
Appendix A. 

A considerable, although purely technical, simplification in the closure scheme outlined 
above based on Eqs. ([5]), (0) occurs when also the attractive potential w(r) in Eq. (0) is 
given by a Yukawa function, i.e. when one has 



w(r) 



exp[— z(r — 1)] 



(8) 



z being the inverse range of the potential. From Eq. 
Eq. (|6|) becomes 



it is then immediately seen that 



g(r) = 



exp[-zi(r- 1)] exp[ 
c(r) = iti h K 2 



Mr ~ 1)] 



r < 1 



r > 1 



(9) 



where K 2 and z 2 are the quantities referred to as K and z in Eq. ([]), and K\, z\ are 
known function of the density. It is now possible to take advantage of the fact that for 
the Ornstein-Zernike equation supplemented by the closure (|9|) extensive analytical results 
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have been determined |24|, [TJJ. If both K\ and K 2 are given, as in the LOGA/ORPA, 



this enables one to solve Eq. @j altogether 26,27]. More generally, irrespective of the 



form of Ki and K 2 , a prescription can be found to determine the reduced compressibility 
Xred as a function of the density p and the internal energy per unit volume u, which can 
be used in Eq. ([J) to obtain a closed PDE. A similar procedure for the same potential 
considered here was adopted in a previous work [|T5]| , where however the hard-sphere 
contribution to the direct correlation function cns( r ) outside the core was not taken into 
account, so that c(r) was given by a simple one- Yukawa tail. This further simplifies the 
theory, but implies that the description of the hard-sphere fluid coincides with that of the 
PY approximation, which as is well known is not very satisfactory at high density. This 
defect becomes more and more severe as the range of the attractive interaction decreases, 
and can considerably affect the phase diagram predicted by the theory, unless some more- 
or-less ad hoc procedure is adopted to correct the hard-sphere thermodynamics. In order 
to incorporate a better treatment of the hard-sphere fluid into the theory one can turn 
to the two- Yukawa form for c(r) of Eq. (^), whose use in the consistency condition ([5]) 
we are now going to illustrate in some detail. In the following we will exploit the results 
determined in Refs. p3|, |24 , [TJJ], which will be respectively referred to as I, II, III. Let us 



introduce the packing fraction £ = np/Q and the quantity 

(10) 




Xred 

which is the square root of the quantity referred to as A in I, II, III. Eq. (|5|) becomes 



(i - (T- \du ) f \dB ) p r \dp*, , 

To obtain a PDE for u we need to express / as a function of p and u in Eq. ((TTJ) . From 
Eq. (11.14) it is found that / can be written as 



, = (*i + Vg(72 ~7i) _ z\ - zl 7172(72 ~ 7i) 

4 [{zi/z 2 ) 72 - {z 2 /z 1 ) Tl ] Zl z 2 [(Vzz) 72 - (Z3A1) 7i] 

where we have set 

The quantities 71 and 72 are given by Eq. (II. 5) 



71 = (14) 



u. 







W 

The ratios U\/Uq and W\/Wq depend on the integrals 



72 = 2-y/q--^. (15) 



/+00 
drr exp[— Zi(r — 1)] g(r) (x = 1,2). (16) 
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^From Eq. (1.35) it is found in fact 



Wi _ 4 + 2z 2 - 4 ~ 1 
W 2(2 + z 2 ) a 2 I 2 - 1 



;i7) 



and the corresponding relation with W\/Wq replaced by Ui/Uq and the index 2 changed 
to 1. The quantities n and Uj depend only on Z{ and are given by Eq. (1.34): 



T; 



1 

2z~ 
1 

2^ 



Zi + 2 
zf + 2z 



+ exp(-Zi 
4 



4 + 2* 



+ exp( 



(19) 



with i — 1, 2. From the expression of the potential (|S|) it is readily seen that I 2 is directly 
related to the internal energy per unit volume u given by Eq. (01): 



(20) 



Eqs. (0), 



allow then to express r y 2 explicitly as a function of p and u: 
A + 2z 2 - z\ 2r 2 u + p 



72 



2{2 + z 2 ) 2a 2 u + p 



(21) 



We now need 71 as a function of p and u. This is less straightforward than for 72, since 
the integral I\ does not have any direct thermodynamic meaning, the exponential in 1\ 
being related to the tail of the direct correlation function of the hard-sphere gas. We have 
then to make use of some further results determined in I — III. ^From Eq. (1.36) it is found 
that the amplitudes Ki, K 2 of the Yukawa functions in the closure @ can be expressed 
in terms of the above introduced quantities U , U\, W , W\. One has 



2{z x + 2) Vf 
K i — Trr^> u o 



Ui 
U 



n 2 



«1 



(22) 



where a\ is given by Eq. (1.37): 



a 1 



(4 + 2z x - zpn 
2(2 + z x )a x 



(23) 



and the corresponding equations with the index 1 replaced by 2 and Uq, U\ replaced by 
Wo, W\. Let us now introduce the quantities x, y given by 



x 



471 

z 2 



(24) 
(25) 
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^From Eq. (III. 30) one has 



U 
W 



X) 



-2 s (Vq-yf 



(26) 
(27) 



where p and s must satisfy Eq. (11.39) (in the notation of II one has x = u q i/u q0 , y 
w ql /w q0 , p = u q0 , s = w q0 ): 

P + s+ _ 9 As o (y-x) 2 



y 2 



±Z 2 



X 



p + s - 



4p 



z l 



(y 



X 



±z 2 - U 2 



{21 



Eq. (]28f) is readily solved for p and s to give 



P 



{Az 2 2 (y - x) 2 - 16y 2 {y - x) 2 - {z\ - z 2 2 ) [z\ - z\ + 4(y 2 - x 2 )] } , (29) 



64 (y — x) 

and the expression for s is obtained by exchanging zi, z 2 and x, y in the r.h.s. of Eq. 
If Eqs. (0), ([2§, (|26|), (0) are used in Eq. (|22|) we finally obtain 



AQ-y/q-aJiy/q 



X 



-{zl-zl) zl-z 2 2+ A{y 2 



{±4(y 



X 



I6y 2 (y - : 
384 iz\ 



(zi + 2) 



zl)a\ 



Ki(y 



x) 



(30) 



and a similar equation obtained by exchanging the indices 1 and 2 and the quantities x, y. 
We recall that in Eq. ( |3"0"D K±, Z\, a±, and a± are known functions of the density p which 
refer to the hard-sphere system. For given values of p and u, Eqs. fl2~l~|), ( p5| ) allow one 
to determine y. Eq. (|3(]) can then be solved numerically with respect to x to obtain 7! 
via Eq. (0). This solves the problem of determining 7! in terms of p and u. The partial 
derivative (df/du) p that appears in Eq. (|Tl|) can then be determined as 



'df 
du 



'df 
dji 



'<hi 

du 



'df_ 
dl2 



'dj2 
du 



(31) 

must be determined 



where {d^ 2 /du) p is calculated explicitly by Eq. (f21|), while {d^i/du) p 
as the derivative of the function implicitly defined by Eq. (|30"D . If we write Eq. ( |30"1) as 
F(x,y,p) = 0, it is found straightforwardly that Eq. ( |TT| ) takes the form 



B(p,u 



du 



C(p,u 



d 2 u 



d(3 ~ vr, ~',9p 2 ' 

where the functions B(p,u) and C(p,u) are given by the following expressions: 



B(p,u) 



V d l2 



;i-£) 2 du 



df dF dx df dF dy 
$72 dx dji d^i dy dj2 



d F dx 



(32) 

(33) 
(34) 



All the partial derivatives in Eqs. (p3|) , (|34| ) are calculated at constant p and can be deter- 
mined by Eqs. (^), fl2lD , (H), (PB|) , (^0|). The resulting expressions are then evaluated as 
a function of p and w via the procedure described above. The same procedure also allows 
one to determine the reduced compressibility as 1/xred = / 2 /(l ~ £) 2 once / nas been 
obtained from Eq. (|12|). The PDE fl32|) is a non-linear diffusion equation that must be 
integrated numerically. To prevent the occurrence of any numerical instability, especially 
in the critical and sub-critical region, we have adopted an implicit finite-differences algo- 



rithm [28 1 tailored to equations that, although globally non-linear, depend on the partial 
derivatives of the unknown function in a linear fashion like Eq. ( P^) . The integration with 
respect to f3 starts at (3 = and goes down to lower and lower temperatures. Before each 
integration step Eq. (|30| ) is solved numerically and the coefficients B(p,u), C(p,u) are 
determined. The density p ranges in a finite interval (0,po), whose high-density bound- 
ary has been typically set at po — 1. The initial condition can be determined by taking 
into account that at (3 = the radial distribution function coincides with that of the 
hard-sphere gas. From Eqs. (^) and (|8[) one has then 

r+co 

u(p, j3 = 0) = — 27cp / drr exp[— z 2 (r — l)]g H s( r ) f° r every p, (35) 

Jo 

where gm( r ) is obtained in the present scheme by the closure (0). For such a closure, 
as shown in Appendix A, both Uq and U\ in Eq. can be determined analytically as 
a function of p, thus providing 7i(p) at (3 = 0. This allows one to obtain u in Eq. fl3~5|) 
analytically as well: in fact, one can solve Eq. (|l|) for 72 as a function of 71, /, and p, 
where / is readily obtained by using the Carnahan-Starling expression of Xred in Eq. ([P0|). 
Once 72 is known, Eq. ([H]) is solved with respect to u. It must be noted that solving 
Eq. fll2|) for 72 gives two branches, so attention must be paid in order to single out the 
branch that actually corresponds to the physical solution. We also need two boundary 
conditions at p = and p = p . ^From Eq. one has immediately 

u(p = 0, (3) — for every (3 . (36) 

At high density we instead make use of the so-called high-temperature approximation 
(HTA), according to which the excess Helmholtz free energy per unit volume is determined 
via Eq. ([35|) for every temperature. In the fluid region of the phase diagram this of course 
is not exact unless (3 = 0, but it becomes more and more accurate as the density of the 
system is increased f29j, so we expect that for a given sweep along the /3-axis the results 
will not differ appreciably from what would be obtained using an hypothetical exact 
boundary condition, provided the boundary p is located at sufficiently high density. We 
used the HTA at p = p for the reduced compressibility. This yields via Eq. (|5|) the 
boundary condition 

fyp(f>o>P) = -q^(P0'P = °) for every /3. (37) 

We have checked that the output of the numerical integration of Eq. fl32|) is quite in- 
sensitive to the specific choice of the high-density boundary condition. Moreover, for 
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Po — 1 moving the boundary condition to higher densities also leaves the results un- 
affected. Eq. ( |l0|) shows that to be physically meaningful, the quantity / has to be 



non-negative. On the other hand, below the critical temperature the solution of Eq. 
does not satisfy this condition along the whole density interval (0,po), but only outside a 
certain temperature-dependent region (p s i((3), p s2 (/?))• For p = p s i((3) or p = p S 2(P) the 
quantity / vanishes, and consequently the compressibility diverges. As (3 changes, p s i(P) 
and p S 2{(3) give then respectively the low- and the high- density branch of the spinodal 
curve predicted by the theory. The fact that / becomes negative for p s i(P) < P < Psiift) 
not only implies that the theory behaves unphysically in this interval, but it also gives rise 
to an analytical instability which would make the numerical integration of the PDE ( p2|) 
impossible, if one tried to determine the solution over the whole interval (0, po) even below 
the critical temperature. Therefore, the region bounded by the spinodal has been excluded 



from the integration of Eq. (j32|) . Specifically, as soon as it is found that / changes sign, 
so that for a certain density p one has f(p,/3) < 0, the integration is restricted to the 
interval (0, p — Ap) or (p + Ap, p ) respectively for p < p c or p > p c , where Ap is the 
spacing of the density grid. Within the precision of the numerical discretization, one has 
Psi = p — Ap (or p s2 = p + Ap) and the further boundary conditions 

u{p si ,(3) = u s (p si ) i = l,2, (3>(3 C , (38) 

where (3 C is the critical inverse temperature and ug{p) is the value of the internal energy 
per unit volume when the compressibility at density p diverges. This can be determined 
by setting / = in Eq. flT2]) and solving for 71 as a function of p and 72. If Eqs. ( ^4|) 
and fl25p are substituted into Eq. (|30| ), an equation for 72 is obtained that allows one to 
determine the value of 72 when 1/Xrcd = for a certain p. Solving Eq. ( pTj ) with respect 
to u then yields Us(p)- 

Once the internal energy per unit volume u has been determined from Eq. (|32"D, the 
pressure P and the chemical potential p are obtained by integration with respect to (3 
via the relations d(f3P)/df3 = — u + pdu/dp, d(/3p)/d(3 = du/dp. Thanks to the self- 
consistency of the theory, this route to the thermodynamics is equivalent to integrating 
the inverse compressibility with respect to p, but it does not require one to circumvent the 
forbidden region in order to reach the high-density branch of the subcritical isotherms. 



3 Results 



The numerical integration of the PDE ( p2|) with the initial condition ( pq ) and the bound- 
ary conditions (|36D-(|38|) has been performed on a density grid with Ap = lCT 3 -iCr 4 . 
At the beginning of the integration the temperature step A/3 was usually set at A/3 = 
2 x 10~ 5 -10~ 5 . As the temperature approaches its critical value, Af3 can be further de- 
creased if one wishes to get very close to the critical point, and then gradually expanded 
back. The integration was usually carried down to (3 ~ 2A(3 C . The inverse range pa- 
rameter of the attractive tail in Eq. (|8|) has been set at z — 1.8. For this value of z 

I 271 predictions have already been 
1 shows the SCOZA results for the compressibility factor 



several simulations [3D, 31, [3^] and theoretical 
reported in the literature. Fig. 
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Z = P/(pkBT) along two different isotherms, corresponding to T = 2 and T = 1.5, to- 



gether with the MC simulation results by Henderson and coworkers |30|. The agreement 



is very good both at low and high density. The compressibility factors are also reported 
in Tab. 1, together with those obtained by the LOGA/ORPA via the internal energy 



route p7| , which is the one that gives the best agreement with the simulation results. 
It can be seen that for the non-critical states reported here the SCOZA and the energy 
route of the LOGA/ORPA are very close to each other. In Tab. 2 the predictions for the 
chemical potential and the reduced compressibility are compared to the data from the 
MC simulations performed in this work. The internal energy per particle is reported in 
Tab. 3, where again the LOGA/ORPA results are also shown. The critical point predicted 
by the theory has been located by the vanishing of the inverse compressibility 1/Xred- No 
extrapolation procedure to 1/Xred = is necessary, since the algorithm adopted here al- 
lows one to get as close as desired to the critical singularity. As mentioned in Sec. 2, 
below the critical temperature T c the theory yields a spinodal curve. The coexistence 
curve must be determined by a Maxwell construction, i.e. by imposing the equilibrium 
conditions p(p g ,T) = p,(pi,T), P(p g ,T) = P(p h T) for the densities p g , pi of the gas and 
liquid phases at coexistence at a temperature T. In comparing our results for the critical 
point and the coexistence curve with the available simulation data, we found that the 
two simulations for the phase diagram of the system under study already reported in the 
literature |3T , |32|j do not agree very well with each other. We then performed a new set 
of simulations using the MC-FSS method summarized in Appendix B. The SCOZA and 
the simulation results for the critical point are compared in Tab. 4, which also shows the 
predictions of other theories ||. It can be seen that the agreement between the SCOZA 
and the present simulation is remarkably good: the error in the critical density and tem- 
perature is respectively slightly more and slightly less than 0.6%. The SCOZA and the 
simulation coexistence curve in the temperature-density plane are compared in Fig. 2. A 
similar comparison in the temperature-internal energy and in the temperature-chemical 
potential plane is shown respectively in Fig. 3 and in Fig 4. In every case the SCOZA 
agrees very well with the simulation. It can be also appreciated that in the SCOZA the 
coexistence curve extends up to the critical point, while, as already observed in the In- 
troduction, this is not always the case with other theories. In Tab. 4 and in Figs. 2-4 
we have also reported the predictions of the simpler version of the SCOZA mentioned in 
Sec. 2, in which the direct correlation function outside the repulsive core is given by just 
one Yukawa tail, and the thermodynamics of the hard-sphere gas is described in the PY 
approximation. It can be observed that even for not very short-ranged interaction the 
treatment of the repulsive contribution considerably affects the phase diagram predicted 
by the theory, the two- Yukawa SCOZA sensibly improving over the one- Yukawa version. 
The behavior of the SCOZA in the critical region has been studied both analytically ]10 



and numerically [|C2 |. This investigation has shown that above the critical temperature the 
SCOZA yields the same critical exponents as the mean spherical approximation (MSA), 
i.e. 7 = 2, 8 — 5, a — — 1, where the usual notation for the critical exponents has been 
used. On the other hand, on the coexistence curve the critical exponents are neither 
spherical nor classical, and one finds 7' = 7/5, a' = —1/10, (3 = 7/20 (here of course (3 
is the critical exponent that gives the curvature of the coexistence curve near the critical 
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point). Although these results were determined in the case of a nearest-neighbor attrac- 
tive lattice gas, we expect them to hold also in the continuum case. Fig. 5 shows the 
reduced compressibility of the HCYF for T > T c and p = p c as a function of the reduced 
temperature t = (T — T c )/T c on a log-log plot. Also shown is the correspondent effective 
exponent 7 e g, defined as the local slope of the plot. It can be seen that 7 e fr eventually 
saturates at 7 = 2, thus signaling the onset of a MSA-like power- law behavior, but the 
asymptotic regime can be detected only at very small reduced temperature (t ~ 10~ 6 ). 
This is the same scenario previously found in the nearest-neighbor lattice gas. For the 
HCYF, the crossover is controlled as expected by the inverse range parameter z. It has 
been verified that as the interaction becomes longer and longer ranged, the asymptotic 



regime is further pushed to smaller and smaller values of the reduced temperature t [33 



4 Conclusions 

We have studied the thermodynamics and the phase diagram of the HCYF using both 
the SCOZA and MC simulations supplemented by a finite-size scaling analysis. The 
comparison between theory and simulation results shows that the SCOZA yields both 
very good overall thermodynamics and a remarkably accurate coexistence curve up to 
the critical point. The version of the SCOZA considered here takes into account the 
hard-sphere contribution to the direct correlation function outside the repulsive core, and 
sensibly improves over the simpler one- Yukawa version, in which the hard-sphere gas is 
described as in the PY approximation. On the other hand, as stated in Sec. 2, here 
(as well as in the simpler version just mentioned) consistency has been enforced between 
the internal energy and the compressibility route, but not between the virial route and 
either of the above. We think that the further development of making the theory fully 
self-consistent by taking also the virial route into account is worth pursuing, since we 
anticipate that the present version of the SCOZA will yield liquid-state pressures from 
the virial theorem that are not significantly better than those obtained using the virial 
theorem with the LOGA/ORPA g(r). We defer a detailed examination of this issue to a 
later study. In this respect it is worth mentioning an investigation of the HCYF along the 
lines considered here [Q, where some results for the critical parameters were reported 
taking into account all the three routes to the thermodynamics although, as explicitly 
pointed out by the authors, the SCOZA equations were studied in an approximate fashion, 
and no attempt to determine the phase diagram was made. 

Although dealing with a Yukawa potential entails certain analytical simplifications in 
implementing the SCOZA, such an approach can be applied to any kind of tail potential. 
It should also be pointed out that the idea of using the requirement of self-consistency 
to get a closed theory of thermodynamics and correlations is pertinent not only to the 
realm of simple fluids or lattice gases, but has also proven to be a powerful tool in the 
study of a system of spins with continuous symmetry |35[], and of a site-diluted |36| or 



random-field []37fl Ising model. 
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A Waisman parameterization of Chs(V 



In this Appendix we recall the procedure that allows one to determine analytically the 
amplitude K\ and the inverse range Z\ of the direct correlation function chs(^) of the hard- 
sphere gas in the Waisman parameterization fl7|). The relevant equations are reported in 
Ref. PH , which will be referred here as IV. Both z\ and K\ are conveniently expressed 
in terms of two quantities Vo, V\ which are formally analogous to Uo, U\ and Wo, W\ 
introduced in Eqs. (|14D , (|T5|). ^From Eq. (IV.2.32 a) it is found that z\ is given by 



q-p 



(V + f - q)f + J(V + P - q)V q 



(39) 



where / and q are defined in Eqs. fllCf), (|T3|). The expression of the amplitude Ki is the 
same as in Eqs. (|22|) with U , U\ replaced by V , V\\ 



2(z 1 +2) V? 



Vn 



(40) 



where ^ is the packing fraction £ = 7rp/6 and a>i is a function of Z\ given by Eq. (|23|). The 
ratio Vi/Vq can be expressed as a function of V by Eqs. (IV. 2. 24) and (IV.2.26). One has 



V 



2Vo^/q 



(Vo + f 2 -q)(V + f 



(41) 



To obtain the explicit expressions of Z\ and Ki as functions of the density, one must then 
feed into Eqs. (|3^)-(|4"T|) the expression of Vo- This depends on the contact value of the 
radial distribution function y = g{r = l + ) via Eq. (IV.2.32 b): 



Vo = QZy -f 2 + l. 



(42) 



For a hard-sphere gas yo can be determined from the equation of state via the virial 
equation: 

— = l + 4£y . (43) 
P 

The requirement that both the virial and the compressibility route to the thermodynam- 
ics must give the Carnahan-Starling equation of state is then satisfied if the Carnahan- 
Starling pressure and compressibility are substituted respectively in Eq. (^) and Eq. (|T0|) . 
Eqs. ( |42| ) and ([41]) then yield V and Vi/Vq as a function of density. From Eqs. ( p9|) 
and ( f40"D we finally get Z\ and K x . 
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B Simulation details 



The principal aspects of the simulation and finite-size scaling techniques employed in this 
work have previously been detailed elsewhere in the context of a similar study of the 
Lennard- Jones fluid. Accordingly we confine our description to the barest essentials and 
refer the reader to reference |T3 for a fuller account of our methods. 

The Monte-Carlo simulations were performed using a Metropolis algorithm within 



the grand canonical ensemble ||38|| . The MC scheme comprises only particle transfer 
(insertion and deletion) steps, leaving particle moves to be performed implicitly as a 
result of repeated transfers. To simplify identification of particle interactions a linked-list 
scheme was employed. This involves partitioning the periodic simulation space of volume 
L 3 into m 3 cubic cells, each of side the cutoff r c . This strategy ensures that interactions 
emanating from particles in a given cell extend at most to particles in the 26 neighbouring 
cells. 

In our Yukawa system the potential was cutoff at a radius r c = 3.0a, and a correction 
term was applied to the internal energy to compensate for the trunction. System sizes 
having m = 3, 4, 5, 6 and 7 were studied, corresponding (at coexistence) to average particle 
numbers of approximately 230,540,1050,1750 and 2900 respectively. For the m = 3,4 
and 5 system sizes, equilibration periods of 10 5 Monte Carlo transfer attempts per cell 
(MCS) were utilised, while for the m = 6 and m = 7 system sizes up to 2 x 10 6 MCS 
were employed. Sampling frequencies ranged from 20 MCS for the m = 3 system to 150 
MCS for the m = 7 system. The total length of the production runs was also dependent 
upon the system size. For the m = 3 system size, 1 x 10 7 MCS were employed, while for 
the m = 7 system, runs of up to 6 x 10 7 MCS were necessary. 

In the course of the simulations, the observables recorded were the particle number 
density p = N/V and the energy density u = E/V. The joint distribution p L (p,u) was 
accumulated in the form of a histogram. In accordance with convention, we express p 
and u in reduced units: p* = pa 3 ,u* = ua 3 . To allow us to explore efficiently the phase 
space of the model, we employed the histogram reweighting technique ||39|| . This method 
allows histogram accumulated at one set of model parameters to be reweighted to provide 
estimates appropriate to another set of not-too-distant model parameters. 

To facilitate study of the subcritical coexistence region, the multicanonical preweight- 
ing technique |3(J was employed. This technique allows one to circumvent the problems 
of metastability and nonergodicity that would otherwise arise from the large free energy 
barrier separating the coexisting phases. Details of this technique and its implementation 



in the fluid context are described in reference [JIG . 

The critical point parameters were estimated using finite-size scaling technique as de- 
scribed in fl6| . This involves matching the distribution function of the ordering operator 
to the independently known universal critical point form appropriate for the Ising univer- 
sality class. The ordering operator is defined as A4 oc (p* + su*), where s is a non-universal 
"field mixing" parameter, which is finite in the absence of particle-hole symmetry, and 
which is chosen to ensure that p(A4) is symmetric in M.. The estimate of the apparent 
critical temperature obtained by this matching procedure is, however, subject to errors 
associated with corrections to finite-size scaling. To deal with this, we extrapolate to the 
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thermodynamic limit using the known scaling properties of the corrections, which are 
expected to diminish (for sufficiently large system sizes) like L' 6 ' v @, where 9 is the 
correction to scaling exponent and v is the correlation length exponent. The extrapola- 
tion has been performed using a least squares fit to the data for the four largest system 
sizes. The results of the extrapolation are shown in figure Al, from which we estimate 
T c = 1.212(2). The associated estimate for the critical density is p* = 0.312(2). 
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TABLE 1 



rp* 


P* 


MC T 


SCOZA 


LOGA/ORPA^ 


OO 


0.4 


2.52 


2.518 


2.518 


oo 


0.6 


4.22 


4.283 


4.283 


OO 


0.8 


7.65 


7.750 


7.750 


2.0 


0.4 


1.08 


1.120 


1.118 


2.0 


0.6 


2.04 


1.977 


1.974 


2.0 


0.8 


4.27 


4.433 


4.432 


1.5 


0.4 


0.69 


0.667 


0.663 


1.5 


0.6 


1.21 


1.220 


1.214 


1.5 


0.8 


3.31 


3.333 


3.330 



Compressibility factor PV/NksT for the hard-sphere + Yukawa fluid (z = 1.8). Density 
and temperature are in reduced units p* = per 3 , T* = k B T/e, where a is the hard-sphere 
diameter and e is the strength of the attractive potential, t: Monte Carlo data from 
|30[. o: LOGA/ORPA-energy route results from Ref. |27[. 



Ref. 
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TABLE 2 





P* 


MC T 


SCOZA 


MC T 


SCOZA 


oo 


0.4 


1.736(2) 


1.7316 


0.1958(2) 


0.19744 


oo 


0.6 


4.833(2) 


4.8147 


0.0848(5) 


0.08721 


2.0 


0.4 


-0.936(2) 


-0.9396 


0.4992(8) 


0.50439 


2.0 


0.6 


0.515(2) 


0.5003 


0.1594(5) 


0.15976 


1.5 


0.4 


-1.823(2) 


-1.8258 


0.968(3) 


1.0150 


1.5 


0.6 


-0.905(2) 


-0.9294 


0.2217(5) 


0.22147 



Chemical potential \i and reduced compressibility x rc d of the hard-sphere Yukawa fluid 
(z = 1.8). Density and temperature are in reduced units, f: Monte Carlo simulation 
performed in this work. The numbers in brackets give the error in the last figure. 
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TABLE 3 





P* 


MC T 




SCOZA 


LOGA/ORPA° 


oo 


0.4 


-2.495 


-2.516(2) 


-2.517 


-2.517 


oo 


0.6 


-3.975 


-4.002(2) 


-4.002 


-4.002 


oo 


0.8 


-5.573 




-5.611 


-5.611 


2.0 


0.4 


-2.583 


-2.595(2) 


-2.583 


-2.574 


2.0 


0.6 


-4.030 


-4.036(2) 


-4.030 


-4.026 


2.0 


0.8 


-5.622 




-5.620 


-5.618 


1.5 


0.4 


-2.622 


-2.640(2) 


-2.623 


-2.602 


1.5 


0.6 


-4.051 


-4.053(2) 


-4.043 


-4.036 


1.5 


0.8 


-5.630 




-5.623 


-5.621 


1.0 


0.6 


-4.073 




-4.097 


-4.065 


1.0 


0.8 


-5.635 




-5.631 


-5.628 



Internal energy per particle of the hard- sphere Yukawa fluid. All quantities are in reduced 
units, f: Monte Carlo simulation of Ref. |3(J. p. Monte Carlo simulation performed in 
this work. The number in brackets give the error in the last figure, o: LOGA/ORPA. 
The entries for T* = oo are from this work, the rest are from Ref. |26| . 
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TABLE 4 





MC° 


MC° 


MC f 


SCOZAV Yuk 


SCOZA 


HMSA* 


MHNC* 


P*c 


0.294 


0.313 


0.312(2) 


0.308 


0.314 


0.36 


0.28 


1 C 


1.192 


1.178 


1.212(2) 


1.201 


1.219 


1.25 


1.21 



Critical density and temperature (in reduced units) for the hard-sphere Yukawa fluid, 
o: MC simulation of Ref. [3T|. O: MC simulation of Ref. |32| f : MC simulation performed 
in this work. p. SCOZA with 1-Yukawa c(r) (see text). •: from Ref. IJ 



19 



FIGURE CAPTIONS 



Fig. 1 Compressibility factor Z = PKpksT) of the hard-sphere Yukawa fluid (z = 1.8) 
as a function of the reduced density p* along two isotherms at reduced temperature 
T* = 2 (upper curve) and T* = 1.5 (lower curve). Full curve: SCOZA. Squares: 
MC simulation results f30fl . 

Fig. 2 Coexistence curve of the hard-sphere Yukawa fluid (z = 1.8) in the density- 
temperature plane. Density and temperature are in reduced units. Full curve: 
SCOZA. Dashed curve: SCOZA with a one- Yukawa direct correlation function c(r) 
(see text). Squares: MC results (this work). 

Fig. 3 Coexistence curve of the hard-sphere Yukawa fluid in the internal energy-temperature 
plane. E*/N is the internal energy per particle in reduced units. Notation as in 
Fig. 2. 

Fig. 4 Coexistence curve of the hard-sphere Yukawa fluid in the temperature-chemical 
potential plane. All quantities are in reduced units. Notation as in Fig. 2. 

Fig. 5 Log-log plot of the reduced compressibility Xrcd of the hard-sphere Yukawa fluid 
(z = 1.8) on the critical isochore as a function of the reduced temperature t = 
(T — T c )/T c according to the SCOZA (a) and effective exponent 7^, defined as 

7cff = -d(logXred)Mlogt) (b). 

Fig. Al The apparent critical temperature, (as defined by the matching condition de- 
scribed in the text), plotted function of L^ 1 )/", with 9 = 0.54 and v = 0.629. 
The extrapolation of the least squares fit to infinite volume yields the estimate 
T* = 1.212(2). 
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